clear
clear matrix
program drop _all

capture log close

cd D:\Projects\TakingtheStreets\Replication			/*Provide appropriate directory information to run the file on your machine*/

log using Logs\GrauvogelLichtVonSoest2016Replication_5_24_16.log, append



set more off
set matsize 1000
set scheme s1mono


*PURPOSE:	This file replicates the primary count analyses reported in Grauvogel,
*			Licht, and von Soest, 2016, "Sanctions and Signals". The program was 
*			written for Stata 13. 
*INSTRUCTIONS:		For smooth operation, be sure to accurately 
*					designate a directory in the <cd> command above. Create  
*					subfolders with the following exact titles in the directory 
*					location:  Logs, Tables, Graphs, and Data. Save the file named
*					GrauvogelLichtVonSoestReplication.dta in the Data folder. The
*					simulation used to produce Fig 1 and Fig A3 may take quite a
*					long time to complete, especialy if your machine is procesisng
*					other tasks. For this reason, the portion of the do-file with
*					the simulation and graphing commands has been commented out. 
*					If you wish to run the simulation and produce the simulated
*					sampling distribution of the quantities of interest, simply
*					delete the /* */ bracketing the end of the do-file.	 
*AUTHOR:	Amanda A. Licht, Binghamton University
*MACHINE:	LICHTBINGOFFICE, Lenovo ThinkStation running StataMP 13
*DATE:		5/24/2016	
*MODIFIED:		
*FILE:		GrauvogelLichtVonSoest2016Replication.do
*LOG:		GrauvogelLichtVonSoest2016Replication_MO_DY_YR.log
*DATA:		GrauvogelLichtVonSoestReplication.dta
*SOURCES:	Protest and Response, Regan 2012
*			Mass Mobilization, phase 2, Clark 2014
*			UCDP/PRIO Armed Conflict v4.0, Themner and Wallensteen 2013
*			COW Formal Alliances v4.1, Gibler 2009
*			Trade Data v3.0, Barbieri and Keshk 2012
*			IGO memberships v2.3, Pevehouse, Nordstrom and Warnke 2004
*			Geddes Autocratic Regimes, Geddes, Wright and Frantz 2013
*			PolityIV Jaggers and Gurr 2012
*			IAEP, Regan and Clark 2013
*			Freedom of the Press, Freedom House 2013
*			CIRI, Cingranelli and Richards 2012
*			Political Terror Scale, Gibney and Wood 2012
*OUTPUT:
*	Data:	GL&vS_CountprobsIVs				- set of mean and modal values of 
*											control variables for use in simulation
*											of quantities of interest (QIs)
*			GL&vS_PredictedProbs&Counts_l0	- simulated sampling distribution of QIs
*											 from Model 1, ZINB, with lag set to 0, 
*											N=100,000
*			GL&vS_PredictedProbs&Counts_l1	- as above, but lag = 1
*			GL&vS_PredictedProbs&Counts_l2	- as above, but lag = 2
*			GL&vS_PredictedProbs&Counts_l5	- as above, but lag = 5
*			GL&vS_PredictedProbs&Counts_CIs	- median, 2.5th and 97.5th percentile of 
*											each QI at each value of Lagged DV
*	Graphs:
*	Tables:	Table1



**************************************
*PRELIMINARIES:
**************************************

*Open up replication data:
use Data/GrauvogelLichtVonSoestReplication.dta
	des
	
*declare the data to be time-series cross-sectional by state and month:
	xtset targetstate yearmonth
	
	
**************************************
*MAIN ANALYSIS:
**************************************	
	
*The following code will run four count models and save the 
*results to a single table with these columns:
*	1) zinb w/lagged DV
*	2) poisson w/fixed effects
*	3) poisson w/fixed effects and a suggested time dummies
*	4) poisson w/fixed effects and lagged DV
*The table will not exactly reflect that reported in the paper,
*largely due to formatting and style choices (e.g. excluding 
*uninteresting parameters, which are reported in the online  
*appendix) but also because the automatic output of the ZINB 
*model will include both equations. Table 1 in the paper only 
*reports the count equation, to emphasize the comparability of
*results across the four models. 


*Model 1, ZINB:																
#delimit ;
zinb MonthlyEventsThatCount 
		threat_thisobs 
		imposition_thisobs 
		HumanRightsSanction 
		NarrowSanction_bin 
		multisender_bin 
		threaton 
		sanctionon
		Loverallglobalizationindex
		politynorm 
		politynorm_sq
		physint
		FOTP_Free 
		FOTP_PFree 	
		Llpop
		LAtrocity
		Increase_stateviolence 
		Increase_stateaccom 
		Decrease_stateviolence 
		Decrease_stateaccom 
		L1MonthlyEventsThatCount, 
	inflate(HumanRightsSanction 
		NarrowSanction_bin 
		multisender_bin
		threaton 
		sanctionon
		Loverallglobalizationindex
		politynorm 
		politynorm_sq
		physint
		FOTP_Free 
		FOTP_PFree
		Llpop
		Llrealgdp
		L1MonthlyEventsThatCount) 
vce(cluster targetstate);
#delimit cr

*Run lincom commands for the implicit interactions, 
*and store the combined coefficients and standard 
*errors in locals to be included in addtext() 
*subcommand of the outreg2 code:
	lincom threaton + HumanRightsSanction	
		local cc1 = round(r(estimate), .001)
		local se1 = round(r(se), .001)
		
	lincom threaton + NarrowSanction_bin
		local cc2 = round(r(estimate), .001)
		local se2 = round(r(se), .001)

	lincom threaton + multisender_bin
		local cc3 = round(r(estimate), .001)
		local se3 = round(r(se), .001)

	lincom sanctionon + HumanRightsSanction	
		local cc4 = round(r(estimate), .001)
		local se4 = round(r(se), .001)
		
	lincom sanctionon + NarrowSanction_bin
		local cc5 = round(r(estimate), .001)
		local se5 = round(r(se), .001)

	lincom sanctionon + multisender_bin
		local cc6 = round(r(estimate), .001)
		local se6 = round(r(se), .001)

*Export table to text file:	
#delimit ;
outreg2 using "Tables\Table1.txt", text  replace
		stats(coef se) aster(coef) sym(***,**,*) cti(ZINB)
		e(N N_zero N_clust ll chi2)
		addtext("Other Threats + HR, `cc1'(`se1')" 
			"Other Threats + Narrow, `cc2'(`se2')"
			"Other Threats + Mult, `cc3'(`se3')"
			"Other Imposed + HR, `cc4'(`se4')"
			"Other Imposed + Narrow, `cc5'(`se5')"
			"Other Imposed + Multi, `cc6'(`se6')");
#delimit cr		

*Save parameter matrices for simulation below:
	mat def b = e(b)
	mat def V = e(V)

	mat list b	

*Model 2, Fixed Effects Poisson:
#delimit ;
xtpoisson MonthlyEventsThatCount 
		threat_thisobs 
		imposition_thisobs 
		HumanRightsSanction 
		NarrowSanction_bin 
		multisender_bin 
		threaton 
		sanctionon
		Loverallglobalizationindex
		politynorm 
		politynorm_sq
		physint
		FOTP_Free 
		FOTP_PFree 	
		Llpop
		LAtrocity
		Increase_stateviolence 
		Increase_stateaccom 
		Decrease_stateviolence 
		Decrease_stateaccom 
		, fe vce(robust);
#delimit cr	
lincom threaton + HumanRightsSanction	
	local cc1 = round(r(estimate), .001)
	local se1 = round(r(se), .001)
	
lincom threaton + NarrowSanction_bin
	local cc2 = round(r(estimate), .001)
	local se2 = round(r(se), .001)

lincom threaton + multisender_bin
	local cc3 = round(r(estimate), .001)
	local se3 = round(r(se), .001)

lincom sanctionon + HumanRightsSanction	
	local cc4 = round(r(estimate), .001)
	local se4 = round(r(se), .001)
	
lincom sanctionon + NarrowSanction_bin
	local cc5 = round(r(estimate), .001)
	local se5 = round(r(se), .001)

lincom sanctionon + multisender_bin
	local cc6 = round(r(estimate), .001)
	local se6 = round(r(se), .001)
	
#delimit ;
outreg2 using "Tables\Table1.txt", text  append
		stats(coef se) aster(coef) sym(***,**,*) cti(FE Poisson)
		e(N N_g ll chi2)
		addtext("Other Threats + HR, `cc1'(`se1')" 
			"Other Threats + Narrow, `cc2'(`se2')"
			"Other Threats + Mult, `cc3'(`se3')"
			"Other Imposed + HR, `cc4'(`se4')"
			"Other Imposed + Narrow, `cc5'(`se5')"
			"Other Imposed + Multi, `cc6'(`se6')");
#delimit cr		



*Model 3, Fixed Effects Poisson w/Time Dummies:
#delimit ;
xtpoisson MonthlyEventsThatCount 
		threat_thisobs 
		imposition_thisobs 
		HumanRightsSanction 
		NarrowSanction_bin 
		multisender_bin 
		threaton 
		sanctionon
		Loverallglobalizationindex
		politynorm 
		politynorm_sq
		physint
		FOTP_Free 
		FOTP_PFree 	
		Llpop
		LAtrocity
		Increase_stateviolence 
		Increase_stateaccom 
		Decrease_stateviolence 
		Decrease_stateaccom 
		 _1992 
		 _1996 
		 _2000 
		 _2001 
		 _2004, fe vce(robust);
#delimit cr	
lincom threaton + HumanRightsSanction	
	local cc1 = round(r(estimate), .001)
	local se1 = round(r(se), .001)
	
lincom threaton + NarrowSanction_bin
	local cc2 = round(r(estimate), .001)
	local se2 = round(r(se), .001)

lincom threaton + multisender_bin
	local cc3 = round(r(estimate), .001)
	local se3 = round(r(se), .001)

lincom sanctionon + HumanRightsSanction	
	local cc4 = round(r(estimate), .001)
	local se4 = round(r(se), .001)
	
lincom sanctionon + NarrowSanction_bin
	local cc5 = round(r(estimate), .001)
	local se5 = round(r(se), .001)

lincom sanctionon + multisender_bin
	local cc6 = round(r(estimate), .001)
	local se6 = round(r(se), .001)
	
#delimit ;
outreg2 using "Tables\Table1.txt", text  append
		stats(coef se) aster(coef) sym(***,**,*) cti(FE Poisson w/Time Dummies)
		e(N N_g ll chi2)
		addtext("Other Threats + HR, `cc1'(`se1')" 
			"Other Threats + Narrow, `cc2'(`se2')"
			"Other Threats + Mult, `cc3'(`se3')"
			"Other Imposed + HR, `cc4'(`se4')"
			"Other Imposed + Narrow, `cc5'(`se5')"
			"Other Imposed + Multi, `cc6'(`se6')");
#delimit cr		


*Model 4, Fixed Effects Poisson w/Lagged DV:
#delimit ;
xtpoisson MonthlyEventsThatCount 
		threat_thisobs 
		imposition_thisobs 
		HumanRightsSanction 
		NarrowSanction_bin 
		multisender_bin 
		threaton 
		sanctionon
		Loverallglobalizationindex
		politynorm 
		politynorm_sq
		physint
		FOTP_Free 
		FOTP_PFree 	
		Llpop
		LAtrocity
		Increase_stateviolence 
		Increase_stateaccom 
		Decrease_stateviolence 
		Decrease_stateaccom 
		L1MonthlyEventsThatCount, fe vce(robust);
#delimit cr	
lincom threaton + HumanRightsSanction	
	local cc1 = round(r(estimate), .001)
	local se1 = round(r(se), .001)
	
lincom threaton + NarrowSanction_bin
	local cc2 = round(r(estimate), .001)
	local se2 = round(r(se), .001)

lincom threaton + multisender_bin
	local cc3 = round(r(estimate), .001)
	local se3 = round(r(se), .001)

lincom sanctionon + HumanRightsSanction	
	local cc4 = round(r(estimate), .001)
	local se4 = round(r(se), .001)
	
lincom sanctionon + NarrowSanction_bin
	local cc5 = round(r(estimate), .001)
	local se5 = round(r(se), .001)

lincom sanctionon + multisender_bin
	local cc6 = round(r(estimate), .001)
	local se6 = round(r(se), .001)
	
#delimit ;
outreg2 using "Tables\Table1.txt", text  append
		stats(coef se) aster(coef) sym(***,**,*) cti(FE Poisson w/Lagged DV)
		e(N N_g ll chi2)
		addtext("Other Threats + HR, `cc1'(`se1')" 
			"Other Threats + Narrow, `cc2'(`se2')"
			"Other Threats + Mult, `cc3'(`se3')"
			"Other Imposed + HR, `cc4'(`se4')"
			"Other Imposed + Narrow, `cc5'(`se5')"
			"Other Imposed + Multi, `cc6'(`se6')");
#delimit cr		



******************************************
*INTERPRETATION
******************************************

*To generate substantively meaningful quantities of interest (QIs) from  
*the ZINB model, the following code will construct a simulated sampling 
*distribution of several QIs. For these caluclations, all other variables
*will be held at their observed means or modes. Therefore, the current
*dataset will be collapsed to obtain these values. 


*Use collapse to get mean values of covariates:
#delimit ;
collapse (mean) Loverallglobalizationindex
				politynorm
				physint
				FOTP_Free 
				FOTP_PFree
				Llpop
				Llrealgdp
				LAtrocity
				Increase_stateviolence 
				Increase_stateaccom 
				Decrease_stateviolence 
				Decrease_stateaccom ;

*For binaries, use mean value to round to modal value:
#delimit ;
foreach var in 	FOTP_Free 
				FOTP_PFree 
				LAtrocity 
				Increase_stateviolence 
				Increase_stateaccom 
				Decrease_stateviolence 
				Decrease_stateaccom		{;
				
	replace `var' = round(`var', 1);
										};
*Since these are just always zero, the XB calculation simplifies quite a bit
tab1 	FOTP_Free 
		FOTP_PFree 
		LAtrocity 
		Increase_stateviolence 
		Increase_stateaccom 
		Decrease_stateviolence 
		Decrease_stateaccom;
#delimit cr

*Save the dataset of control variable values:
	save 	Data/GL&vS_CountprobsIVs.dta, replace	

	
*Write probram to calculate QIs:
	
set more off
*program drop countprobs
program define countprobs														/*Open program*/
	syntax, lag(real)															/*Set syntax to control value of lagged DV*/
	
		use Data/GL&vS_CountprobsIVs.dta, clear									/*Open datset of mean/modal values*/
		expand 1000																/*Duplicate 1000 times to increase sample size*/
		
		*draw vector of betas from the parameter matrices of Model 1
		drawnorm double b1-b37, means(b) cov(V)
		
		
		*Get Alpha Inverse:
		gen alphainv = 1/(exp(b37))
		
		
		*Generate PAllzero, Mu (i.e. the exponentiated XB for the count), 
		*predicted probability of a count greater than zero, & predicted  
		*counts across new threats/impositions given no ongoing:
		
		*First, we need to calculate the linear index for the inflate equation:
		#delimit ;
		gen xb_inflate = 	b27*Loverallglobalizationindex+
							b28*politynorm +
							b29*politynorm*politynorm+
							b30*physint+
							b33*Llpop+
							b34*Llrealgdp+
							b35*`lag'+
							b36;						
		#delimit cr
		
		*Probability of all zero is logit:
			gen Pallzero = 1/(1+exp(-xb_inflate))
	
		*Now get linear indices for each of the types of new sanction, with 
		*the indicators on (1) and off (0). NT will indicate "New Threat";
		*NI, "New Imposition"
		foreach new in 0 1	{													/*Open loop over on/off*/
		
			*Calculate linear index:
			#delimit ;															/*XB|New threat*/
			gen xb_NT`new' =b1*`new' +
							b8*Loverallglobalizationindex+
							b9*politynorm +
							b10*politynorm*politynorm+
							b11*physint+
							b14*Llpop+
							b20*`lag'+
							b21;


			
			gen xb_NI`new' =b2*`new' +											/*XB|New imposed*/
							b8*Loverallglobalizationindex+
							b9*politynorm +
							b10*politynorm*politynorm+
							b11*physint+
							b14*Llpop+
							b20*`lag'+
							b21;
			#delimit cr 
			
			*Mu is the exponentiated linear index:
			gen mu_NT`new' = exp(xb_NT`new')
			
			gen mu_NI`new' = exp(xb_NI`new')
			
			*Calculate probability of a zero count given inflate=0 using ZINB distribution:
			gen ProbZ_NT`new'= (1-Pallzero)*(exp(lngamma(0+alphainv)) / ( round(exp(lnfactorial(0)),1) * exp(lngamma(alphainv)) )  )* ((alphainv/(alphainv+mu_NT`new'))^alphainv) * ((mu_NT`new'/(alphainv+mu_NT`new'))^0)
			
			gen ProbZ_NI`new'= (1-Pallzero)*(exp(lngamma(0+alphainv)) / ( round(exp(lnfactorial(0)),1) * exp(lngamma(alphainv)) )  )* ((alphainv/(alphainv+mu_NI`new'))^alphainv) * ((mu_NI`new'/(alphainv+mu_NI`new'))^0)
			
			*Calculate probability of nonzero count as 1-probability(Y=0|inflate=0)
			gen ProbNZ_NT`new' = 1-ProbZ_NT`new'
			
			gen ProbNZ_NI`new' = 1-ProbZ_NI`new'
			
			*Calculate predicted count:
			gen ECount_NT`new' = (1-Pallzero)*mu_NT`new'
			
			gen ECount_NI`new' = (1-Pallzero)*mu_NI`new'
			
			
				*Label variables:
					label var ProbZ_NT`new' "Pr(Y=0|NT=`new')"
					label var ProbZ_NI`new' "Pr(Y=0|NI=`new')"
					label var ProbNZ_NT`new' "Pr(Y>0|NT=`new')"
					label var ProbNZ_NI`new' "Pr(Y>0|NI=`new')"
					label var ECount_NT`new' "E(Y|NT=`new')"
					label var ECount_NI`new' "E(Y|NI=`new')"
							}													/*Close loop over on/off for new sanctions*/
							
			
		*Now, get these same quantities, but over the values of:  
		*Other Threats, Other Impositions, Targeted, Human Rights, 
		*& Multilateral. I'll do this with a bunch of loops.
		*The naming convention for the predicted quantities is
		*		0	-> "off"
		*		1	-> "on"
		*		NT 	-> new threat
		*		NI	-> new imposed
		*		OT	-> other threat 
		*		OI	-> other imposed
		*		HR	-> human rights issue
		*		narr-> targeted sanction
		*		mult-> multisender sanction
		

		foreach new in 0 1	{													/*Open loop over new sanctions on/off*/
		foreach OT in 0 1	{													/*Open loop over other threats on/off*/
		foreach OI in 0 1	{													/*Open loop over other imposed on/off*/
		foreach HR in 0 1	{													/*Open loop over human rights on/off*/
		foreach narr in 0 1	{													/*Open loop over targeted on/off*/
		foreach mult in 0 1	{													/*Open loop over multisender on/off*/
		
			*Calculate Mu as exponentiated linear index of the count: 
			#delimit ;
			gen mu_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = exp(b1*`new' +
															b2*0 +
															b3*`HR' +
															b4*`narr' +
															b5*`mult' +
															b6*`OT' +
															b7*`OI' +
															b8*Loverallglobalizationindex+
															b9*politynorm +
															b10*politynorm*politynorm+
															b11*physint+
															b14*Llpop+
															b20*`lag'+
															b21);


						
			gen mu_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = exp(b1*0 +
															b2*`new' +
															b3*`HR' +
															b4*`narr' +
															b5*`mult' +
															b6*`OT' +
															b7*`OI' +
															b8*Loverallglobalizationindex+
															b9*politynorm +
															b10*politynorm*politynorm+
															b11*physint+
															b14*Llpop+
															b20*`lag'+
															b21);												
			
			*Calculate linear index of inflate equation:
			#delimit ;
			gen xb_inflate`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = 	b22*`HR' +
																b23*`narr' +
																b24*`mult'+
																b25*`OT' +
																b26*`OI'+
																b27*Loverallglobalizationindex+
																b28*politynorm +
																b29*politynorm*politynorm+
																b30*physint+
																b33*Llpop+
																b34*Llrealgdp+
																b35*`lag'+
																b36;	
				
			*Gen probability of al zero:	
			#delimit ;
			gen Pallzero`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = 	1/(1+exp(-xb_inflate`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'));												
			
			*generate probability of a zero count given inflate=0
			#delimit ;
			gen ProbZ_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'= 
				(1-Pallzero`new'_`OT'`OI'_H`HR'_N`narr'_M`mult')*(exp(lngamma(0+alphainv)) / ( round(exp(lnfactorial(0)),1) * exp(lngamma(alphainv)) )  
				)* ((alphainv/(alphainv+mu_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'))^alphainv) * (
				(mu_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'/(alphainv+mu_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'))^0);
	
			gen ProbZ_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'= 
				(1-Pallzero`new'_`OT'`OI'_H`HR'_N`narr'_M`mult')*(exp(lngamma(0+alphainv)) / ( round(exp(lnfactorial(0)),1) * exp(lngamma(alphainv)) )  
				)* ((alphainv/(alphainv+mu_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'))^alphainv) * (
				(mu_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'/(alphainv+mu_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'))^0);
			#delimit cr
			
			*Calculate probability of nonzero count as 1-probability(Y=0|inflate=0)
			gen ProbNZ_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = 1-ProbZ_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'
			
			gen ProbNZ_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = 1-ProbZ_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'
			
			*Calculate predicted count:
			gen ECount_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = (1-Pallzero)*mu_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'
			
			gen ECount_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' = (1-Pallzero)*mu_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'												

				label var ProbZ_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' "P(Y=0|NT=`new',OT=`OT',OI=`OI',HR=`HR',narr=`narr',multi=`multi')"
				label var ProbZ_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' "P(Y=0|NI=`new',OT=`OT',OI=`OI',HR=`HR',narr=`narr',multi=`multi')"
				label var ProbNZ_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' "P(Y>0|NT=`new',OT=`OT',OI=`OI',HR=`HR',narr=`narr',multi=`multi')"
				label var ProbNZ_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' "P(Y>0|NI=`new',OT=`OT',OI=`OI',HR=`HR',narr=`narr',multi=`multi')"
				label var ECount_NT`new'_`OT'`OI'_H`HR'_N`narr'_M`mult'	"E(Y|NT=`new',OT=`OT',OI=`OI',HR=`HR',narr=`narr',multi=`multi')"
				label var ECount_NI`new'_`OT'`OI'_H`HR'_N`narr'_M`mult' "E(Y|NI=`new',OT=`OT',OI=`OI',HR=`HR',narr=`narr',multi=`multi')"
							}
							}
							}
							}
							}
							}													/*Close all loops*/
							
							
			*Calculate Differences in Probability of a nonzero count across conditions:
			
			*1) new threat/imposition on vs. off
				
				gen DNZ_NT10	= ProbNZ_NT1 - ProbNZ_NT0
				gen DNZ_NI10	= ProbNZ_NI1 - ProbNZ_NI0
				
				gen rDNZ_NT10	= DNZ_NT10/ProbNZ_NT0
				gen rDNZ_NI10	= DNZ_NI10/ProbNZ_NI0
				
			*2) ongoing threats/impositions vs no episodes by type:
	
				gen DNZ_OT10_HR = ProbNZ_NT0_10_H1_N0_M0 - ProbNZ_NT0_00_H0_N0_M0
				gen DNZ_OT10_N 	= ProbNZ_NT0_10_H0_N1_M0 - ProbNZ_NT0_00_H0_N0_M0
				gen DNZ_OT10_M 	= ProbNZ_NT0_10_H0_N0_M1 - ProbNZ_NT0_00_H0_N0_M0
				
				gen DNZ_OI10_HR = ProbNZ_NT0_01_H1_N0_M0 - ProbNZ_NT0_00_H0_N0_M0
				gen DNZ_OI10_N = ProbNZ_NT0_01_H0_N1_M0 - ProbNZ_NT0_00_H0_N0_M0
				gen DNZ_OI10_M = ProbNZ_NT0_01_H0_N0_M1 - ProbNZ_NT0_00_H0_N0_M0
				
					label var DNZ_OT10_HR "Diff in Pr(Y>0) OT HR - no episode"
					label var DNZ_OT10_N "Diff in Pr(Y>0) OT narr - no episode"
					label var DNZ_OT10_M "Diff in Pr(Y>0) OT multi - no episode"
	
					label var DNZ_OI10_HR "Diff in Pr(Y>0) OI HR - no episode"
					label var DNZ_OI10_N "Diff in Pr(Y>0) OI narr - no episode"
					label var DNZ_OI10_M "Diff in Pr(Y>0) OI multi - no episode"
				
				gen rDNZ_OT10_HR = DNZ_OT10_HR/ProbNZ_NT0_00_H0_N0_M0
				gen rDNZ_OT10_N  = DNZ_OT10_N/ProbNZ_NT0_00_H0_N0_M0
				gen rDNZ_OT10_M  = DNZ_OT10_M/ProbNZ_NT0_00_H0_N0_M0
				
				gen rDNZ_OI10_HR = DNZ_OI10_HR/ProbNZ_NT0_00_H0_N0_M0
				gen rDNZ_OI10_N  = DNZ_OI10_N/ProbNZ_NT0_00_H0_N0_M0
				gen rDNZ_OI10_M  = DNZ_OI10_M/ProbNZ_NT0_00_H0_N0_M0
				
					label var rDNZ_OT10_HR "RelativeDiff in Pr(Y>0) OT HR - no episode"
					label var rDNZ_OT10_N "RelativeDiff in Pr(Y>0) OT narr - no episode"
					label var rDNZ_OT10_M "RelativeDiff in Pr(Y>0) OT multi - no episode"
	
					label var rDNZ_OI10_HR "RelativeDiff in Pr(Y>0) OI HR - no episode"
					label var rDNZ_OI10_N "RelativeDiff in Pr(Y>0) OI narr - no episode"
					label var rDNZ_OI10_M "RelativeDiff in Pr(Y>0) OI multi - no episode"

			*3) ongoing threats/impositions  of type vs. generic type:
				gen DNZ_OT_HR10 = ProbNZ_NT0_10_H1_N0_M0 - ProbNZ_NT0_10_H0_N0_M0
				gen DNZ_OT_N10 	= ProbNZ_NT0_10_H0_N1_M0 - ProbNZ_NT0_10_H0_N0_M0
				gen DNZ_OT_M10 	= ProbNZ_NT0_10_H0_N0_M1 - ProbNZ_NT0_10_H0_N0_M0
				
				
				gen DNZ_OI_HR10 = ProbNZ_NT0_01_H1_N0_M0 - ProbNZ_NT0_01_H0_N0_M0
				gen DNZ_OI_N10 	= ProbNZ_NT0_01_H0_N1_M0 - ProbNZ_NT0_01_H0_N0_M0
				gen DNZ_OI_M10 	= ProbNZ_NT0_01_H0_N0_M1 - ProbNZ_NT0_01_H0_N0_M0

					label var DNZ_OT_HR10 "Diff in Pr(Y>0) OT HR - episode"
					label var DNZ_OT_N10 "Diff in Pr(Y>0) OT narr - episode"
					label var DNZ_OT_M10 "Diff in Pr(Y>0) OT multi - episode"
	
					label var DNZ_OI_HR10 "Diff in Pr(Y>0) OI HR - episode"
					label var DNZ_OI_N10 "Diff in Pr(Y>0) OI narr - episode"
					label var DNZ_OI_M10 "Diff in Pr(Y>0) OI multi - episode"				
				
				gen rDNZ_OT_HR10 = DNZ_OT_HR10/ProbNZ_NT0_10_H0_N0_M0
				gen rDNZ_OT_N10  = DNZ_OT_N10/ProbNZ_NT0_10_H0_N0_M0
				gen rDNZ_OT_M10  = DNZ_OT_M10/ProbNZ_NT0_10_H0_N0_M0
				
				
				gen rDNZ_OI_HR10 = DNZ_OI_HR10/ProbNZ_NT0_01_H0_N0_M0
				gen rDNZ_OI_N10  = DNZ_OI_N10/ProbNZ_NT0_01_H0_N0_M0
				gen rDNZ_OI_M10  = DNZ_OI_M10/ProbNZ_NT0_01_H0_N0_M0	
				
					label var rDNZ_OT_HR10 "RelativeDiff in Pr(Y>0) OT HR - episode"
					label var rDNZ_OT_N10 "RelativeDiff in Pr(Y>0) OT narr - episode"
					label var rDNZ_OT_M10 "RelativeDiff in Pr(Y>0) OT multi - episode"
	
					label var rDNZ_OI_HR10 "RelativeDiff in Pr(Y>0) OI HR - episode"
					label var rDNZ_OI_N10 "RelativeDiff in Pr(Y>0) OI narr - episode"
					label var rDNZ_OI_M10 "RelativeDiff in Pr(Y>0) OI multi - episode"
					
		*4) new threats given ongoing threats/impositions of type vs no ongoing
			
				gen DNZ_NT1_OT10_HR = ProbNZ_NT1_10_H1_N0_M0 - ProbNZ_NT1_00_H0_N0_M0
				gen DNZ_NT1_OT10_N 	= ProbNZ_NT1_10_H0_N1_M0 - ProbNZ_NT1_00_H0_N0_M0
				gen DNZ_NT1_OT10_M 	= ProbNZ_NT1_10_H0_N0_M1 - ProbNZ_NT1_00_H0_N0_M0
				

				gen DNZ_NT1_OI10_HR = ProbNZ_NT1_01_H1_N0_M0 - ProbNZ_NT1_00_H0_N0_M0
				gen DNZ_NT1_OI10_N 	= ProbNZ_NT1_01_H0_N1_M0 - ProbNZ_NT1_00_H0_N0_M0
				gen DNZ_NT1_OI10_M 	= ProbNZ_NT1_01_H0_N0_M1 - ProbNZ_NT1_00_H0_N0_M0

					label var DNZ_NT1_OT10_HR "Diff in Pr(Y>0) NT&OT HR-no episode"
					label var DNZ_NT1_OT10_N "Diff in Pr(Y>0) NT&OT narr-no episode"
					label var DNZ_NT1_OT10_M "Diff in Pr(Y>0) NT&OT multi-no episode"
	
					label var DNZ_NT1_OI10_HR "Diff in Pr(Y>0) NT&OI HR-no episode"
					label var DNZ_NT1_OI10_N "Diff in Pr(Y>0) NT&OI narr-no episode"
					label var DNZ_NT1_OI10_M "Diff in Pr(Y>0) NT&OI multi-no episode"	
					
				gen rDNZ_NT1_OT10_HR = DNZ_NT1_OT10_HR/ProbNZ_NT1_00_H0_N0_M0
				gen rDNZ_NT1_OT10_N  = DNZ_NT1_OT10_HR/ProbNZ_NT1_00_H0_N0_M0
				gen rDNZ_NT1_OT10_M  = DNZ_NT1_OT10_HR/ProbNZ_NT1_00_H0_N0_M0
				

				gen rDNZ_NT1_OI10_HR = DNZ_NT1_OI10_HR/ProbNZ_NT1_00_H0_N0_M0
				gen rDNZ_NT1_OI10_N  = DNZ_NT1_OI10_N/ProbNZ_NT1_00_H0_N0_M0
				gen rDNZ_NT1_OI10_M  = DNZ_NT1_OI10_M/ProbNZ_NT1_00_H0_N0_M0

					label var rDNZ_NT1_OT10_HR "RelativeDiff in Pr(Y>0) NT&OT HR-no episode"
					label var rDNZ_NT1_OT10_N "RelativeDiff in Pr(Y>0) NT&OT narr-no episode"
					label var rDNZ_NT1_OT10_M "RelativeDiff in Pr(Y>0) NT&OT multi-no episode"
	
					label var rDNZ_NT1_OI10_HR "RelativeDiff in Pr(Y>0) NT&OI HR-no episode"
					label var rDNZ_NT1_OI10_N "RelativeDiff in Pr(Y>0) NT&OI narr-no episode"
					label var rDNZ_NT1_OI10_M "RelativeDiff in Pr(Y>0) NT&OI multi-no episode"		
					
		*5) new threats given ongoing threats/impositions of type vs generic ongoing type:
				gen DNZ_NT1_OT_HR10 = ProbNZ_NT1_10_H1_N0_M0 - ProbNZ_NT1_10_H0_N0_M0
				gen DNZ_NT1_OT_N10 = ProbNZ_NT1_10_H0_N1_M0 - ProbNZ_NT1_10_H0_N0_M0
				gen DNZ_NT1_OT_M10 = ProbNZ_NT1_10_H0_N0_M1 - ProbNZ_NT1_10_H0_N0_M0
				
				
				gen DNZ_NT1_OI_HR10 = ProbNZ_NT1_01_H1_N0_M0 - ProbNZ_NT1_01_H0_N0_M0
				gen DNZ_NT1_OI_N10 = ProbNZ_NT1_01_H0_N1_M0 - ProbNZ_NT1_01_H0_N0_M0
				gen DNZ_NT1_OI_M10 = ProbNZ_NT1_01_H0_N0_M1 - ProbNZ_NT1_01_H0_N0_M0		
				
					label var DNZ_NT1_OT_HR10 "Diff in Pr(Y>0) NT&OT HR-episode"
					label var DNZ_NT1_OT_N10 "Diff in Pr(Y>0) NT&OT narr-episode"
					label var DNZ_NT1_OT_M10 "Diff in Pr(Y>0) NT&OT multi-episode"
	
					label var DNZ_NT1_OI_HR10 "Diff in Pr(Y>0) NT&OI HR-episode"
					label var DNZ_NT1_OI_N10 "Diff in Pr(Y>0) NT&OI narr-episode"
					label var DNZ_NT1_OI_M10 "Diff in Pr(Y>0) NT&OI multi-episode"	
					
				gen rDNZ_NT1_OT_HR10 = DNZ_NT1_OT_HR10/ProbNZ_NT1_10_H0_N0_M0
				gen rDNZ_NT1_OT_N10 = DNZ_NT1_OT_N10/ProbNZ_NT1_10_H0_N0_M0
				gen rDNZ_NT1_OT_M10 = DNZ_NT1_OT_M10/ProbNZ_NT1_10_H0_N0_M0
				
				
				gen rDNZ_NT1_OI_HR10 = DNZ_NT1_OI_HR10/ProbNZ_NT1_01_H0_N0_M0
				gen rDNZ_NT1_OI_N10 = DNZ_NT1_OI_N10/ProbNZ_NT1_01_H0_N0_M0
				gen rDNZ_NT1_OI_M10 = DNZ_NT1_OI_M10/ProbNZ_NT1_01_H0_N0_M0		
				
					label var rDNZ_NT1_OT_HR10 "RelativeDiff in Pr(Y>0) NT&OT HR-episode"
					label var rDNZ_NT1_OT_N10 "RelativeDiff in Pr(Y>0) NT&OT narr-episode"
					label var rDNZ_NT1_OT_M10 "RelativeDiff in Pr(Y>0) NT&OT multi-episode"
	
					label var rDNZ_NT1_OI_HR10 "RelativeDiff in Pr(Y>0) NT&OI HR-episode"
					label var rDNZ_NT1_OI_N10 "RelativeDiff in Pr(Y>0) NT&OI narr-episode"
					label var rDNZ_NT1_OI_M10 "RelativeDiff in Pr(Y>0) NT&OI multi-episode"	

			*Trim to just the variables that we need:		
			keep Prob* Pallzero* D* rD* ECount*
			
			*Make a variable to hold the value of the lag in this iteration:
				gen LagCount = `lag'
				
			*Save results into a dataset for QIs at this value of lag:
			append using Data/GL&vS_PredictedProbs&Counts_l`lag'.dta
			save Data/GL&vS_PredictedProbs&Counts_l`lag'.dta, replace
			
end		/*close the program*/		


*make datasets for saving the results inside the program:
clear
foreach lag in 0 1 2 5 	{
save Data/GL&vS_PredictedProbs&Counts_l`lag'.dta, replace emptyok
						}

*Make a place to store the confidence intervals, that will have a variable
*for the lagged count values so that values can be merged instead of appending.
*This allows for a nice, neat dataset without excess missing observations:
clear
set obs 4	/*I'll do 10 values of the lag, since this is going to be an intense loop*/
gen LagCount =_n-1
	replace LagCount = 5 in 4
tab LagCount	
save  Data/GL&vS_PredictedProbs&Counts_CIs.dta, replace emptyok	

*Execute the program with 100 repetitions at every value of lagged count.
*This may take a considerable amount of time, especially if your computer is
*doing other work at the same time. Uncomment this block of code to run the 
*simulations and graph the results.


set more off
disp "Time Date:  $S_TIME		$S_DATE"
foreach lag in  0 1 2 5 	{
	simulate, reps(100): countprobs, lag(`lag')
							}
disp "Time Date:  $S_TIME		$S_DATE"

											
*Use statsby to grab the 2.5th, 50th, and 97.5th percentiles of each QI in each
*dataset:
set more off
foreach lag in 0 1 2 5 	{														/*loop over datasets*/
use Data/GL&vS_PredictedProbs&Counts_l`lag'.dta, clear							/*open simulated sampling distribution*/
#delimit ;
foreach var in	DNZ_NT10														/*loop over all the QIs*/
				DNZ_NI10
				DNZ_OT10_HR
				DNZ_OT10_N
				DNZ_OT10_M
				DNZ_OI10_HR
				DNZ_OI10_N
				DNZ_OI10_M
				DNZ_OT_HR10
				DNZ_OT_N10
				DNZ_OT_M10
				DNZ_OI_HR10
				DNZ_OI_N10
				DNZ_OI_M10
				DNZ_NT1_OT10_HR
				DNZ_NT1_OT10_N
				DNZ_NT1_OT10_M
				DNZ_NT1_OT10_HR
				DNZ_NT1_OT10_N
				DNZ_NT1_OT10_M
				DNZ_NT1_OI10_HR
				DNZ_NT1_OI10_N
				DNZ_NT1_OI10_M
				DNZ_NT1_OT_HR10
				DNZ_NT1_OT_N10
				DNZ_NT1_OT_M10
				DNZ_NT1_OI_HR10
				DNZ_NT1_OI_N10
				DNZ_NT1_OI_M10
				rDNZ_NT10
				rDNZ_NI10
				rDNZ_OT10_HR
				rDNZ_OT10_N
				rDNZ_OT10_M
				rDNZ_OI10_HR
				rDNZ_OI10_N
				rDNZ_OI10_M
				rDNZ_OT_HR10
				rDNZ_OT_N10
				rDNZ_OT_M10
				rDNZ_OI_HR10
				rDNZ_OI_N10
				rDNZ_OI_M10
				rDNZ_NT1_OT10_HR
				rDNZ_NT1_OT10_N
				rDNZ_NT1_OT10_M
				rDNZ_NT1_OT10_HR
				rDNZ_NT1_OT10_N
				rDNZ_NT1_OT10_M
				rDNZ_NT1_OI10_HR
				rDNZ_NT1_OI10_N
				rDNZ_NT1_OI10_M
				rDNZ_NT1_OT_HR10
				rDNZ_NT1_OT_N10
				rDNZ_NT1_OT_M10
				rDNZ_NT1_OI_HR10
				rDNZ_NT1_OI_N10
				rDNZ_NT1_OI_M10	
				ProbZ_NT1  
				ProbZ_NT0
				ProbZ_NI1  
				ProbZ_NI0
				ProbZ_NT0_00_H0_N0_M0
				ProbZ_NT0_10_H0_N0_M0
				ProbZ_NT0_01_H0_N0_M0
				ProbZ_NT1_00_H0_N0_M0
				ProbZ_NT1_01_H0_N0_M0
				ProbZ_NT1_10_H0_N0_M0
				ProbZ_NT0_10_H1_N0_M0
				ProbZ_NT0_10_H0_N1_M0
				ProbZ_NT0_10_H0_N0_M1
				ProbZ_NT0_01_H1_N0_M0
				ProbZ_NT0_01_H0_N1_M0
				ProbZ_NT0_01_H0_N0_M1
				ProbZ_NT1_10_H1_N0_M0
				ProbZ_NT1_10_H0_N1_M0
				ProbZ_NT1_10_H0_N0_M1
				ProbZ_NT1_01_H1_N0_M0
				ProbZ_NT1_01_H0_N1_M0
				ProbZ_NT1_01_H0_N0_M1
				ProbNZ_NT1  
				ProbNZ_NT0
				ProbNZ_NI1  
				ProbNZ_NI0
				ProbNZ_NT0_00_H0_N0_M0
				ProbNZ_NT0_10_H0_N0_M0
				ProbNZ_NT0_01_H0_N0_M0
				ProbNZ_NT1_00_H0_N0_M0
				ProbNZ_NT1_01_H0_N0_M0
				ProbNZ_NT1_10_H0_N0_M0
				ProbNZ_NT0_10_H1_N0_M0
				ProbNZ_NT0_10_H0_N1_M0
				ProbNZ_NT0_10_H0_N0_M1
				ProbNZ_NT0_01_H1_N0_M0
				ProbNZ_NT0_01_H0_N1_M0
				ProbNZ_NT0_01_H0_N0_M1
				ProbNZ_NT1_10_H1_N0_M0
				ProbNZ_NT1_10_H0_N1_M0
				ProbNZ_NT1_10_H0_N0_M1
				ProbNZ_NT1_01_H1_N0_M0
				ProbNZ_NT1_01_H0_N1_M0
				ProbNZ_NT1_01_H0_N0_M1
				ECount_NT1  
				ECount_NT0
				ECount_NI1  
				ECount_NI0
				ECount_NT0_00_H0_N0_M0
				ECount_NT0_10_H0_N0_M0
				ECount_NT0_01_H0_N0_M0
				ECount_NT1_00_H0_N0_M0
				ECount_NT1_01_H0_N0_M0
				ECount_NT1_10_H0_N0_M0
				ECount_NT0_10_H1_N0_M0
				ECount_NT0_10_H0_N1_M0
				ECount_NT0_10_H0_N0_M1
				ECount_NT0_01_H1_N0_M0
				ECount_NT0_01_H0_N1_M0
				ECount_NT0_01_H0_N0_M1
				ECount_NT1_10_H1_N0_M0
				ECount_NT1_10_H0_N1_M0
				ECount_NT1_10_H0_N0_M1
				ECount_NT1_01_H1_N0_M0
				ECount_NT1_01_H0_N1_M0
				ECount_NT1_01_H0_N0_M1
				Pallzero0_00_H0_N0_M0
				Pallzero0_10_H0_N0_M0
				Pallzero0_01_H0_N0_M0
				Pallzero1_00_H0_N0_M0
				Pallzero1_01_H0_N0_M0
				Pallzero1_10_H0_N0_M0
				Pallzero0_10_H1_N0_M0
				Pallzero0_10_H0_N1_M0
				Pallzero0_10_H0_N0_M1
				Pallzero0_01_H1_N0_M0
				Pallzero0_01_H0_N1_M0
				Pallzero0_01_H0_N0_M1
				Pallzero1_10_H1_N0_M0
				Pallzero1_10_H0_N1_M0
				Pallzero1_10_H0_N0_M1
				Pallzero1_01_H1_N0_M0
				Pallzero1_01_H0_N1_M0
				Pallzero1_01_H0_N0_M1
				Pallzero	{;

	use Data/GL&vS_PredictedProbs&Counts_l`lag'.dta, clear;						/*open simulated sampling distribution*/

	*Use statsby to grab correct percentiles:
	#delimit ;
	statsby  `var'_lo = r(r1) `var'_med = r(r2) `var'_hi = r(r3), clear by(LagCount): 
			_pctile `var', p(2.5, 50, 97.5);
	#delimit cr
	
	*Merge median and confidence bounds with CIs dataset by value of LagCount
	merge 1:1 LagCount using Data/GL&vS_PredictedProbs&Counts_CIs.dta, nogen
	save Data/GL&vS_PredictedProbs&Counts_CIs.dta, replace
											}									/*close loop over variables*/
															}					/*close loop over datasets*/

															
************************
*Graph the results:
************************

*One thing we'll want to do is make two variables to serve as the X-axis for 
*graphing the differences. These will be different from lagcount, so that 
*the values will be equally spaces. Also, one of them will be slightly larger, 
*so that the results for new impositions can be "jittered" over to the side, 
*increasing readability.

sort LagCount
gen Xlag = _n - 1

	lab def Xfmt 0 "0" 1 "1" 2 "2" 3 "5" 
	lab val Xlag Xfmt
	tab Xlag
	

gen Xlag2 = Xlag + .2
		label var Xlag2 "Xlag + .2 for graphing clarity"
		
*results are also more visually appealing if graphed as percentages rather than decimals:
foreach var in 	rDNZ_NT10_lo rDNZ_NT10_med rDNZ_NT10_hi	rDNZ_NI10_lo rDNZ_NI10_med rDNZ_NI10_hi	{
	gen `var'pct = `var'*100
																								}

*Replicate Figure 1:																								
#delimit ;
twoway 	rspike rDNZ_NT10_lopct rDNZ_NT10_hipct Xlag , sort lc(black)
		||rspike rDNZ_NI10_lopct rDNZ_NI10_hipct Xlag2 , sort lc(gs10)
		||scatter rDNZ_NT10_medpct Xlag , sort msym(square) mc(black) msi(small)
		||scatter rDNZ_NI10_medpct Xlag2 , sort msym(square) mc(gs10) msi(small)
			xlabel(0(1)3, valuelabel) 
			xtitle("Lagged Count of Protests")
			ytitle("% Change in Pr(Y>0)")
			yline(0, lc(gs12) lp(dash))
	title("Fig 1. Relative Differences in Predicted Probability of Nonzero Count of" "Antigovernment Protests by New Sanctions and Previous Levels", size(medium))		
	legend(order(3 "New Threats" 4 "New Impositions") ring(0) pos(7) cols(1) size(small))
	note("NOTE: Spikes provide 95% credible intervals from simulated sampling distribution of 100,000 draws from the parameter matrices"
		"matrices of the ZINB reported as Model 1 in Table 1. All covariates held at mean/modal values for the calculation. Relative"
		"difference is the Pr(Y>0|sanction) - Pr(Y>0|~sanction)/Pr(Y>0|~sanction).", size(vsmall) span)
	name(fig1, replace)
	saving(Graphs/Fig1.gph, replace);
#delimit cr	


*Replicate Figure A3:
#delimit ;
twoway 	rspike ECount_NT1_lo ECount_NT1_hi Xlag if LagCount<=5, sort lc(black)
		||rspike ECount_NI0_lo ECount_NI0_hi Xlag2 if LagCount<=5, sort lc(gs10)
		||scatter ECount_NT1_med Xlag if LagCount<=5, sort msym(circle) mc(black) msi(small)
		||scatter ECount_NI0_med Xlag2 if LagCount<=5, sort msym(square) mc(gs10) msi(small)
			xlabel(0(1)3, valuelabel) 
			xtitle("Lagged Count of Protests")
			ytitle("Predicted Count of Protests")
	title("Threat of Sanction", size(small))		
	legend(order(3 "Yes" 4 "No") title("New Threat?", size(small)) ring(0) pos(11) cols(1) size(small))
	name(EctT, replace);
#delimit cr			

#delimit ;
twoway 	rspike ECount_NI1_lo ECount_NI1_hi Xlag if LagCount<=5, sort lc(black)
		||rspike ECount_NI0_lo ECount_NI0_hi Xlag2 if LagCount<=5, sort lc(gs10)
		||scatter ECount_NI1_med Xlag if LagCount<=5, sort msym(circle) mc(black) msi(small)
		||scatter ECount_NI0_med Xlag2 if LagCount<=5, sort msym(square) mc(gs10) msi(small)
			xlabel(0(1)3, valuelabel) 
			xtitle("Lagged Count of Protests")
			ytitle("")
	title("Imposition of Sanction", size(small))		
	legend(order(3 "Yes" 4 "No") title("New Imposition?", size(small)) ring(0) pos(11) cols(1) size(small))
	name(EctI, replace);
#delimit cr		

#delimit ;
graph combine EctT EctI, ycommon title("Fig A3. Predicted Count of Antigovernment Protests by New Sanction Event" "over previous levels", size(medium))
	note("NOTE:Spikes indicate 95% credible intervals from simulated sampling distribution of 100,000 draws from the variance-covariance" 
		"matrices of the model with controls reported in Table 2. All covariates held to mean or modal values for the calculations.", size(vsmall) span)
	name(figa3, replace)
	saving(Graphs/FigA3.gph, replace);
#delimit cr


*This completes the main analysis. Proceed to GrauvogelLichtVonSoest2016Replication2.do for robustness checks,
*autocorrelation analysis, and bivariate probits.  
log close
